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Abstract 

We address the question concerning the maximum density, /O^axi which nuclei 
(and more generally - nuclear structures) can exist in neutron star interiors. An ab¬ 
solute upper bound to p^ax is obtained using the bulk approximation, in which sur¬ 
face and Coulomb effects are neglected. A very good approximation to p^ax is given 
by the threshold for the instability of a uniform npe plasma with respect to density 
modulations; this threshold is calculated using the Extended Thomas-Fermi approxi¬ 
mation for the Skyrme energy functionals. For recent SLy Skyrme forces, which are 
particularly suitable for the description of very neutron rich nucleon systems, one gets 
Pmax = 0-08 fm“^; at this density protons constitute only 4% of nucleons. 

PACS numbers: 97.60. Jd, 21.65. 95.30. Cq 


1 Introduction 

Neutron drip instability limits the neutron excess of neutron-rich nuclei which can be formed 
in laboratory to d = {N — Z)IA < dn-drip — 0.3. However, this limitation is no longer valid 
for nuclei in the interiors of neutron stars, where the density varies from a few g cm“^ near 
the stellar surface to more than 10^^ g cm“^ near the star center. Above 10^ g cm“^, atomic 
structures are crushed, and electrons form an essentially uniform Fermi gas. Standard 
scenario of neutron star formation in gravitational collapse of a massive stellar core (or of 
a mass accreting white dwarf) predicts nuclear composition corresponding to the state of 
complete thermodynamic equilibrium (minimum of the free energy per nucleon). Neutron 
excess parameter of nuclei, immersed in electron gas, increases with increasing density (i.e., 
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increasing depth below stellar surface). Beta decay of neutron rich nuclei, which would be 
unstable in vacuum, is blocked via Pauli exclusion principle due to the presence of dense, 
degenerate electron gas. Up to the density 4 —6 x 10^^ g cm“^, neutron star matter consists 
of nuclei immersed in dense electron gas. At higher densities, neutrons start to populate 
continuum states, forming a degenerate neutron gas. The presence of an outer degenerate 
neutron gas influences the properties of nuclei by: a) blocking (due to Pauli exclusion 
principle) neutron emission from nuclei, b) exerting pressure on nuclei (i.e. compressing 
them ), and c) modifying (lowering) nuclear surface energy. 

Further increase of density is accompanied by an increase of the fraction of volume 
occupied by nuclei and a simultaneous decrease of proton fraction in neutron star matter 
(and in nuclei). In classical terms, nucleon component of neutron star matter consists 
there of two coexisting nucleon fluids: denser one in the interior of nuclei, and the less 
dense outer neutron gas (at highest densities the less dense nucleon fluid can contain some 
admixture of protons, see Section 2). At some density the two-fluid phase becomes 
unstable with respect to transition into a uniform, electrically neutral npe plasma, composed 
mostly of neutrons, with a few percent admixture of electrons and protons. The density 
is the maximum density, at which nuclei can exist in the neutron star interiors; it 
turns out to be significantly lower than normal nuclear saturation density po = 0.16 fm“^, 
typically ~ (see Section 3). Before p^ax i® reached, the interplay between Coulomb and 
surface effects can lead to the appearance of unusual shapes of ‘nuclei’ (rods, plates, tubes, 
bubbles, ...). However, both the theoretical value of p^^x the actual shape of ‘nuclei’ 
(the more appropriate term would be ‘nuclear structures’) at highest densities, depend on 
the assumed effective nuclear hamiltonian, used in many-body calculations i ii0|]- 
For some effective N-N interactions, neutron star matter contained spherical nuclei down 
to p()(ax (SkM force in Q, all relativistic mean field models in |^). For other models of 
effective N-N forces, spherical nuclei were replaced by a sequence of cylindrical ones (rods), 
flat (plates), cylindrical holes in nuclear matter filled by neutron gas (tubes), and finally 
spherical holes in nuclear matter filled by neutron gas (bubbles), before disappearence of 
nuclear structures at p^ax [0; 1; 0) Hi • 

Let us stress that up to about po the elementary constituents of neutron star matter are 
still the same as those of terrestrial matter: n, p, e (see, e.g., [|^). However, under extremal 
conditions prevailing in neutron star interior, the structures these constituents form can be 
dramatically different from those familiar from terrestrial ones. 

The quantity p^^x has an important astrophysical meaning: it determines the bottom 
boundary of the neutron star crust, an outer envelope, in which nuclei form a periodic 
crystal lattice due to longe range Coulomb interactions. In the case of unusual nuclei, one 
deals rather with two dimensional (rods, tubes) or even one dimensional (plates) liquid 
crystals il 10- The solid crust plays an important role in the evolution and dynamics 
of neutron stars (see, e.g. |13|). 

In view of the fact, that the bottom layers of neutron star crust involve nucleon matter 
with an extremely large neutron excess (5 ~ 0.9), an appropriate choice of effective nuclear 
hamiltonian is of crucial importance. Recently, a new set of the Skyrme-type effective 
N-N interactions SLy (Skyrme Lyon) has been derived, based on an approach which is 
particularly appropriate, as far as the applications to a neutron rich matter are concerned 
[0, Relevant additional experimental items concerning nuclei with large neutron excess. 
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isovector effective masses, constraints of spin stability, and requirement of consistency with 
the realistic UV14+VIII equation of state (EOS) of dense neutron matter of Wiringa et 
ah 0 for pq < p < l.b iva. were included into a fitting procedure for the SLy forces 
parameters. In the present paper we calculate using the SLy models of effective N-N 
interaction, and compare our results with those obtained using older Skyrme-type forces. 
Ski' and SkM*, used frequently in astrophysical applications. The parameters of the SLy 
forces, used in the present calculations, together with those of the SkM* and Ski' models, 
are given in Table 1. The SLy4 is a basic SLy force; the SLy7 model has been obtained 
following the most ambitious fitting procedure, in which spin-gradient terms and center of 
mass correction term were both included in the Skyrme energy functional |Q]. 

At the densities of interest, matter in the interior of a neutron star which is more than 
one year old, is strongly degenerate, and thermal contributions to thermodynamic quantities 
can be safely neglected. In what follows, we will consider the properties of neutron star 
matter in the T = 0 approximation; dense matter will be assumed to be in its ground state 
(it is then called ‘cold catalyzed matter’). 

Table 1 

Parameter values of the Skyrme forces 


force 

SLy4 

SLy7 

SkM* 

SkT 

to (MeV fm^) 

-2488.91 

-2482.41 

-2645.0 

-1057.3 

ti (MeV fm^) 

486.82 

457.97 

410.0 

235.9 

t 2 (MeV fm^) 

-546.39 

-419.85 

-135.0 

-100.00 

h (MeV fm3+3'') 
a 

Xq 

13777.0 

1 

6 

0.834 

13677.0 

1 

6 

0.846 

15595.0 

1 

6 

0.09 

14463.5 

1 

0.2885 

Xl 

-0.344 

-0.511 

0 

0 

X2 

-1.000 

-1.000 

-1.000 

0 

X3 

1.354 

1.391 

0 

0.2257 

Wo (MeV fm^) 

123.0 

126.0 

130.0 

120.0 


In Section 2, we present the calculation of p^^^ in the bulk approximation. Threshold 
for the instability of homogeneous npe matter with respect to density modulations (i.e., 
formation of ‘nuclear structures’) is calculated in Section 3. Properties of very neutron 
rich nuclei in the bottom layers of neutron star crust, and dependence on the effective N-N 
interaction model, are discussed in Section 4. Concluding remarks are presented in Section 
5. 
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2 Bound for the existence of nuclei: the bulk approximation 


Above neutron drip density, Pn-drip) neutron star crust is a two-phase nucleon system, the 
denser nucleon phase (fluid) ‘i’ residing inside nuclei and the less dense ‘o’ one forming a gas 
outside nuclei. Both nucleon fluids are permeated by an essentially homogeneous electron 
gas, which ensures overall charge neutrality. The shape of nuclei, i.e. that of the ‘i-o’ 
interface, results from the balance of the surface and Coulomb terms in the total energy of 
the system. 

The bulk approximation consists in neglecting the Coulomb and surface effects. At given 
mean nucleon density p, nucleons are present in general in both of the coexisting ‘i’ and ‘o’ 
fluids, characterized by the constant densities pni, Ppi, Pno, and above proton drip density, 
/?p_drip, also Ppo- The equilibrium between the ‘i’ and ‘o’ fluids, ensured by the (strong) N-N 
interaction, implies the equality of the chemical potentials of nucleons, and the equality of 
nucleon pressures, 


P-ni — pnoi Ppi — Ppo ; Pni — Pno j (1) 

where the label N indicates the nucleon contribution to the matter pressure, and the condi¬ 
tion on the proton chemical potentials applies only above the proton-drip density, Pp-drip- 
At a given p, Eqs. i) enables determination of pni, Ppi, Pno, and for p > pp_drip, also ppo- 
Let us denote the volume fraction occupied by the ‘i’ fluid by u. At a given mean nucleon 
density p = up\ + {1 — u)po, the total energy density is given by 

E = uEjqi -|- (1 — u)E}^o + Eg , (2) 

where Eg is the electron energy density, Eg{pg). Beta equilibrium within the npe matter 
implies relation between the chemical potentials of neutrons, protons and electrons. 


pn — pp T Pe ) (3) 

while requirement of the overall charge neutrality leads to 

Pe = ttppi -k (1 - u)ppo . (4) 

Eqs. (HJ^^) determine completely the equilibrium of a two-fluid npe matter at a given 
nucleon density p. 

The two-fluid phase is stable at a given p if its energy per nucleon is lower than that in 
the uniform (one-fluid) phase (i.e., that of a uniform npe matter). This stability condition 
breaks down at a density P 2 ^i- Calculations show that approaching p 2 ^i from the lower 
density (two-fluid) side corresponds to u —> 1, i.e. to the ‘i’ (denser) fluid filling the whole 
volume jl^] . Therefore the 2 <-> 1 transition is a continuous one, with no density jump 
(our results for the SLy4 force are shown in Eig. 1). Actually, Coulomb and surface effects 
would increase the energy per nucleon in the two-fluid phase, as compared to the plain bulk 
approximation. Therefore, the real transition into a uniform npe matter occurs at density 
lower than p 2 «i, so that p 2 ^i is thus a robust upper bound to p(^ax- 
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Figure 1. Composition and the two-fluid - one-fluid phase transition in the bulk approxi¬ 
mation, for the SLy4 force. Two-fluid phase to the left, one-fluid phase to the right of the 
dotted line representing the bottom boundary of the crust in the bulk approximation. 
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Table 2 

Neutron drip and proton drip densities, and proton fraction 
and nucleon density at the 2-fluid^ 1-fluid transition, 
calculated in the bulk approximation 


force 

Pn—drip 

(fm“^) 

Pp—drip 

(fm“^) 

7*2^ 1 

(fm“^) 

(%) 

SkM* 

7.55 X 10"^ 

0.0781 

0.0866 

3.13 

SkT 

8.09 X 10-^ 

0.1066 

0.1088 

3.71 

SLy4 

7.34 X 10"'^ 

0.0853 

0.0905 

3.84 

SLy7 

7.21 X lO-"^ 

0.0836 

0.0889 

3.80 


Numerical values of p 2 ^i for four Skyrme forces are shown in Table 2, where we show 
also the proton-drip densities, and the proton fractions at the 2 1 transition point, 

X 2 ^i (bulk instability of two-fluid phase with Skyrme forces SkM and Ski' was studied in 
[^]). . As far as the value of p 2 ^i is concerned, both SLy models yield values close to 
0.09 fm“^, intermediate between those corresponding to the Ski' and SkM* forces, which 
yield 0.11 fm“^ and 0.087 fm“^, respectively. Proton drip takes place in a narrow interval 
of densities close to p 2 ^i- At the 2 1 transition point, protons constitute less than 4% 

of nucleons. 


Instability of a uniform npe matter and existence of nuclear 
structures 


An instability of a spatially uniform state of the npe matter with respect to density mod¬ 
ulations signals the appearence of nuclear structures in the true ground state of the npe 
system. Let us consider first a npe matter at nucleon density p > po, which is sufficiently 
high to exlude any possibility of existence of nuclear structures. By decreasing p, we will 
eventually find the threshold density, below which uniform npe matter is unstable with re¬ 
spect to density modulations. This threshold density turns out a lower bound (and a rather 
good approximation) to p^^^ 


0 


The ground state of uniform npe matter of a given nucleon density p corresponds to the 
minimum of the energy density E{pn, Pp, Pe) = Eq, calculated at a fixed nucleon density 
p = Pp+Pn = p, under the condition of charge neutrality pe = Pp (all densities being assumed 
constant in space). Minimisation implies the beta equilibrium relation between the chemical 
potentials of matter constituents, pn = Pp + Pe- This relation ensures vanishing of the first 
variation of E implied by small perturbations Spj{r) (where j = n,p, e). However, this does 
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not guarantee the stability of the spatially homogeneous state of the npe matter, which 
requires that the second variation of E (quadratic in 5pj) be positive. 

The energy per unit volume, E, of a slightly nonuniform npe matter, can be decomposed 
into nucleon, electron, and Coulomb components. 


E — En + T'e + .E'Coul • 


(5) 


The nucleon contribution to E can be expressed in terms of the energy density 

SNiPn, Pp,^Pn, Vpp), obtained from the Skyrme models of the effective nucleon hamiltonian, 

so that for a perturbed, slightly spatially inhomogeneous state we get the energy functional 

En = ^ J drSN[pn{r), Pp{r),\7pnir),Vpp{r)] . ( 6 ) 


The expression for £]\f has been calculated using the semi-classical Extended Thomas-Fermi 
(ETF) treatment of the kinetic and the spin-orbit energy densities |§]. Assuming that 
the spatial gradients are small, one keeps only the quadratic gradient terms in the ETF 
expressions. This approximation will be justified by the fact that characteristic wavelengths 
of perturbations turn out to be much larger then the internucleon distance (see below). With 
such approximations, the change of the energy (per unit volume) implied by the density 
perturbations can be expressed, keeping only second order terms, as in, 


E-En = I 


dq 

{2ttY 


■Fij{q}6pi{q)Spj{qy 


(7) 


where Fourier representation 


6pj{r) = 


dq 

{2TrY 


6pj{q)e 


iqr 


( 8 ) 


has been used. The hermitian F)j(q) matrix determines the stability of the uniform state 
of equilibrium of the npe matter with respect to the spatially periodic perturbations of the 
densities of wavevector q. Due to the isotropy of the homogeneous equilibrium state of the 
npe matter, Fij depends only on |q| = q. In the case of Skyrme forces, the matrix elements 
Eij can be calculated analytically, as explicit functions of the equilibrium densities and q, 
and are composed of compression (local), curvature (gradient), and Coulomb components 


0- 


The condition for the Eij matrix to be positive definite turns out to be equivalent to 
the requirement that the determinant of the Eij matrix be positive Q. At each density 
p, one has thus to check the condition det[Fjj(( 7 )] > 0. Let us start with some p > pq, at 
which det[F)j(g)] > 0 for any q. By decreasing p, we find eventually a wavenumber Q at 
which stability condition is violated for the first time; this happens at some density p{Q)- 
For p < p{Q) homogeneous state is no longer the true ground state of the npe matter since 
it is unstable with respect to small periodic density modulations. 

In contrast to the bulk approximation, in which constant densities were assumed by 
construction, the general energy functional used in this section allows for a consistent treat¬ 
ment of spatially inhomogeneous states of the npe matter. The instability at p{Q) signals 
a phase transition with a loss of translational symmetry of the npe matter. In principle. 
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this could be a second order phase transition. However, the combination of Coulomb and 
surface effects turns out to be sufficiently strong to destabilize the npe state with an in- 
hnitesimal density modulation, leading to a first-order phase transition into a state with 
finite amplitude density modulations. The real equilibrium phase transition (at constant 
pressure) will thus take place at a npe matter density pi, slightly higher than p{Q). Homo¬ 
geneous npe matter of density pi coexists with a spatially inhomogeneous phase exhibiting 
some nuclear structures, of a density p 2 slightly lower than pi. Therefore, p{Q) is actually 
a lower bound on pi, but in view of the closeness of pi and p 2 , it can be used as a rather 
good approximation of = P2- 


Table 3 

Threshold for the instability of the homogeneous npe matter 
with respect to density modulations 


force 

p{Q) 

(fm“^) 

(fm“^) 

x{Q) 

(%) 


a 

b 

c 

a 

b 

c 

c 

SkM* 

0.0738 

0.0744 

0.0754 

0.277 

0.284 

0.299 

2.79 

SkT 

0.0992 

0.0993 

0.1005 

0.314 

0.319 

0.367 

3.52 

SLy4 

0.0781 

0.0787 

0.0794 

0.262 

0.271 

0.281 

3.57 

SLy7 

0.0773 

0.780 

0.0786 

0.270 

0.280 

0.292 

3.55 


a: TF approximation without spin-orbit term, 
b: ETF without spin-orbit term, 
c: ETF including spin-orbit term. 


Our numerical results for four Skyrme forces are shown in Table 3 ([|l^], the case of 
Ski' force was studied previously, using slightly different approximations, in j^). Column 
“a” was obtained using standard TF approximation to nucleon kinetic energy densities, 
and neglecting spin-orbit contribution to Sn. Adding quadratic gradient terms in the ETF 
expansion of the nucleon energy densities (columns “b”) increased only slightly the values 
of p{Q) (by less than 1%) as compared to the TF values. Notice that both “a” and “b” 
were obtained neglecting the spin-orbit term in Finally, adding the spin-orbit term, 
calculated using the ETF approximation (column “c”), further increased p{Q) by about 









1% with respect to the “b” ones (from 0.8% for SLy7 to 1.5% for Ski'). In the case of the 
critical wavenumber Q, corresponding effects are significantly higher. 

For the ETF approximation to be correct, the value of the characteristic wavelength of 
critical density perturbations, 27r jQ, must be significantly larger than the mean internucleon 
distance. Despite a small proton fraction, 27r jQ ~ 17 — 22 fm is typically four times higher 
than the distance between protons Vpp = (47rp/3)“^/^; for neutrons this ratio is typically 
about eight. 

4 Nuclei in the bottom layers of the crust 

The actual structure of the bottom layer of neutron star crust should be calculated including 
Coulomb and surface effects. One usually considers a limited set of five possible nuclear 
shapes: spheres of i-phase embedded in o-gas, rods of i-phase in o-gas, equidistant plates 
of i-phase with o-gas between them, tubes in the i-matter filled with o-gas, and spherical 
bubbles of the o-gas in i-matter. For a given nuclear shape, the ground state of the matter 
is calculated using the Wigner-Seitz approximation, in which the real system is replaced by 
a set of non-interacting, electrically neutral cells, each cell containing one nucleus. Wigner- 
Seitz cells are spherical, cylindrical or plate-like, depending on the symmetry of nuclear 
structure. The total volume of cells is equal to the volume of the real system. 

By comparing the enthalpies per nucleon {E + P)/p at given P, one hnds the actual 
shape of nuclei at this pressure P. Depending on the assumed effective N-N interaction, 
unusual nuclear shapes (starting with rods) may appear as early as at p ~ 5 — 6 x 10“^ fm“^ 
(the case of FPS force of i, i, 0), or do not appear at all, spherical nuclei being present 
down to p 2 (SkM force case of |^, [^). Elementary considerations within the liquid drop 
model indicate, that at lower densities spherical shape is the only possible one. A possibility 
(or rather a necessity) of appearance of unusual nuclear shapes at highest densities results 
from the fact, that the fraction occupied by i-phase (i.e., nuclear matter), u, increases with 
increasing p, and above some limiting density the sum of the surface and Coulomb energies 
can be reduced by changing nuclear shape from spherical to an unusual one. The main 
point is whether unusual shape will appear before the transition into a homogeneous npe 
matter. 

In general, the structure of the bottom layers of the neutron star crust may be expected 
to be rather sensitive to the behavior of the nuclear surface tension, a, in the relevant 
density interval ^po < p < |po. Our results obtained within the ETF approximation for 
the Skyrme energy functional, shown in Fig. 2 visualize strong dependence of a at large 
neutron excess on the Skyrme force used [notice that all these forces lead to very similar 
values of a{5i = 0), consistent with experimental value extracted from the liquid droplet 
model]. At = 0.8 (which corresponds to the proton fraction Xi = 0.1), characteristic of 
p ~ 0.08 fm“^, nuclear surface tension for SkM* model is only half of that obtained for the 
Ski' one, which in turn is 30% higher than our values for the SLy forces. We should remind, 
however, that even at p ~ |po we have to deal with nuclei which are very far from the most 
neutron rich nuclei available in terrestrial nuclear physics, and therefore we should rely on 
the extrapolation of nuclear models to very high neutron excess. This may be visualized by 
our preliminary results obtained within the compressible liquid drop model for SLy4 force, 
shown in Table 4. 
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Figure 2. Surface tension of the i-o phase interface versus neutron excess parameter (5i = 
1 — 2x\ in the denser (i) phase far from the interface. Typical values of {x\) range from 
0.6 (0.2) at p = 0.05 fm-^ to 0.8 (0.1) at p = 0.08 fm-^. 
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Table 4 


Examples of spherical nuclei in the bottom layers of neutron star crust 
with SLy4 force. Curvature corrections to surface energy are neglected, 
^ceii is the number of nucleons in Wigner-Seitz cell, u is the fraction 
of volume occupied by protons, Z and A are numbers of protons and nucleons 
in nuclei, Rn, Rp are corresponding (equivalent) radii. 


? 

CO 

Pol Pi 

u 

Rp 

(fm) 

Rn Rp 

(fm) 

Z 

A 

^cell 

0.0160 

0.0932 

0.0154 

5.56 

0.83 

28.2 

133.2 

751.3 

0.0425 

0.2819 

0.0551 

6.20 

0.86 

25.1 

154.9 

769.3 

0.0501 

0.3441 

0.0748 

6.42 

0.85 

24.3 

162.5 

744.8 

0.0628 

0.4584 

0.1269 

6.98 

0.80 

23.6 

186.9 

705.6 


5 Conclusion 

Neutron star interiors are expected to contain extremely neutron rich nuclei (or more gen¬ 
erally “nuclear structures”, including those with unusual shapes), with neutron excess far 
beyond the laboratory neutron-drip limit. However, no stable nuclear structure can exist 
above a specific limiting density, which turns out to be significantly smaller than the satu¬ 
ration density of symmetric nuclear matter. For Skyrme forces SLy, which are particularly 
suitable for applications in neutron star matter calculations, this ultimate density is about 
0.08 fm“^, the proton fraction at the limiting density being about 4%. 

The procedure of construction of SLy forces make them suitable also for the calculation 
of equation of state of neutron star matter at supranuclear densities, up to Spo @1- Therefore 
these effective N-N interactions can be used for a unified description of the whole neutron 
star, including its liquid interior. A low value of implies a small mass of neutron star 
crust. Consider a neutron star model of a “canonical mass” 1.4 Mq (measured masses of 
binary radio pulsars are quite close to this canonical value). For SLy7 force, the mass of 
the crust constitutes then only 1.3% of stellar mass, while contribution of the crust to the 
total moment of inertia of neutron star is about 2.8% |]l^]. 

In spite of its small mass, the outer mantle of neutron star, containing neutron rich 
nuclei - neutron star crust - is of paramount importance for star evolution and dynamics. 
It insulates thermally neutron star surface from its liquid interior, and therefore plays 
an essential role in neutron star cooling. It can accumulate mechanical stresses, leading 
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to instabilities crucial for the glitches in pulsar timing. Finally, neutron star crust can 
support, in contrast to the liquid interior, nonaxial deformations (“mountains”) which, 
combined with rapid rotation, could make neutron star a source of continuous gravitational 
radiation. This gives a strong astrophysical motivation for theoretical studies of neutron 
rich nuclei in neutron star interiors. 
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